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0\ • ABSTRACT 

a^ ■ 

We present a new shear estimator for weak lensing observations which properly accounts 
^ ' for the effects of a realistic point spread function (PSF). Images of faint galaxies are subject 

to gravitational shearing followed by smearing with the instrumental and/or atmospheric PSF. 
We construct a 'finite resolution shear operator' which when applied to an observed image has 
the same effect as a gravitational shear applied prior to smearing. This operator allows one 
to calibrate essentially any shear estimator. We then specialize to the case of weighted second 
moment shear estimators. We compute the shear polarizability which gives the response of an 
' individual galaxy's polarization to a gravitational shear. We then compute the response of the 

■ population of galaxies, and thereby construct an optimal weighting scheme for combining shear 

estimates from galaxies of various shapes, luminosities and sizes. We define a figure of merit — 
' an inverse shear variance per unit solid angle — which characterizes the quality of image data 

0^ ■ for shear measurement. The new method is tested with simulated image data. We discuss the 

I correction for anisotropy of the PSF and propose a new technique involving measuring shapes 

i-^ ' from images which have been convolved with a re-circularizing PSF. We draw attention to a 

hitherto ignored noise related bias and show how this can be analyzed and corrected for. The 
O I analysis here draws heavily on the properties of real PSF's and we include as an appendix a 

' brief review, highlighting those aspects which are relevant for weak lensing. 



Subject headings: gravitational lensing - dark matter - clusters of galaxies - large-scale structure 



1. Introduction 

Weak lensing provides a probe of the dark matter distribution on a range of scales from galaxy halos, 
through clusters of galaxies, to large-scale-structure (see e.g. the recent review of [Mellier (1999| ), and 



references therein). In the weak lensing or thin lens approximation, the effect of gravitational lensing on the 
image of a distant object is a mapping of the surface brightness: 

nr.)=fii6^,-^.,)r,) (1) 

where the 2-vector is the angular position on the sky measured relative to the center of the image, / is 
the intrinsic surface brightness that would be seen in the absence of lensing, and the symmetric 'distortion 
tensor' is an integral along the line of sight of the transverse second derivatives of the gravitational 



potential $ (Gimn 1967). In an open cosmology, for instance, the distortion for an object at conformal 



distance can be written as 



, , sinhcij sinhfo^s — Cij) „ „ ^ 

V'/™ = 2 / duj -did^^ (2) 

smhw. 
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(Bar-Kana 1996; Kaiser 1998| ) where di = d/dxi and Xi — Oisinhuj with 0i being a 2-component Cartesian 



vector in the plane of the sky and where the potential is related to the density contrast by V^^ = AnGSp, 
the Laplacian here being take with respect to proper spatial coordinates. Equation can be generalized 
to deal with sources at a range of distances, and with cither accurately known redshifts or partial redshift 
information from broadband colors. 

The distortion will, in general, change both the shapes and sizes, and hence luminosities, of distant 
objects. Any component of the distortion which is coherent over large scales — larger than the typical 
angular separation of background galaxies — is therefore potentially observable as a relative modulation of 
the counts of objects or as a statistical anisotropy of the galaxy shapes, and this allows one to constrain 
the fluctuations in the total density 6p on the corresponding scales. Here we will focus on analysis of 
shape anisotropy, or 'image shear', though the methodology is readily extendible to include the effects of 
magnification. 

While the effect of a shear on the sky surface brightness (|l|) is rather simply stated, no completely 
satisfactory method for estimating the distortion has yet emerged. Perhaps ideally one would attack 
this problem using likelihood; that is, one would ask: what is the probability to observe a given set of 
background galaxy images given that they are drawn from some statistically isotropic unlensed parent 
distribution, as a function of the parameters tpim- Unfortunately, this does not seem to be a particularly 
tractable problem. Also, the problem is further complicated by the finite resolution of real observations, 
and by noise in the images. Instead, what has been done is to adopt some plausible shape statistic — 
typically some kind of central second moment — and then compute how this responds to a gravitational 
shear. We will now review the various different shear estimators that have been proposed, how they relate 
to one another, and what their limitations are. 



1.1. Projection Matrices and the Shear Operator 



As a preliminary, we now introduce some mathematical formalism which will simplify the analysis. 
Symmetric 2x2 tensors like ipim feature prominently in what follows. Such tensors have 3 real degrees of 
freedom. For instance, it is conventional to parameterize the three real degrees of freedom of ipim by the 
triplet comprising the convergence k and the shear 7q,, a = 1, 2 with 



^3 



(3) 



H + ll 72 

A simple way to convert between triplet and tensor components is to use the three constant 2x2 'projection 
matrices': 



Mciij — 



1 
1 



M 



lij 





-1 



1 

1 



(4) 



The symmetrized products of pairs of these matrices [AIaMb] = \{MAimMBmn + J^AnmMBmi) have 
multiplication table 

" Af Ml M2 ' 

[MaMb] = Ml Mo = SabMo + {SaoSbc. + 5aJb^)M^ (5) 
_M2 Mo. 

from which follow the useful identities that the contractions of products and triple products are 

MAlmMBml = 25 AB (6) 
MAlmMBmnMcnl = '2.{5 Bc5 AO + 5 Ac5 BO + 5 AbSco) ■ (7) 
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Any symmetric tensor i/,„ can be written as a linear combination of projection matrices with coefficients 

tA, that is, tira = tAMAlm, and using (|) we have MBlmUm = MBlmMAlmtA = 2SABtA, so tA = jMAlmtlm- 

In this language the convergence and the shear are the three components of the triplet representation of 
the distortion tensor: k = ^M^^imipim — i'o and 7^ = ■^Maim^'im = V'qj = li^. We will adopt the 
convention that upper case Roman indices range over 0, 1, 2 while lower case Greek symbols range over 1, 2 
and that repeated indices are to be summed over. A two component 'polar' ta transforms under rotations as 
ta — > Ra[}{'26)tfj while to transforms as a scalar under rotations. An alternative and widely used formalism 



( Schneider, Ehlers, fc Falco 1992 ) is to regard 71, 72 as the real and imaginary parts of a complex shear, 
but we shall not adopt that approach here. 

It is also convenient to define a 'shear operator' S-y, which generates the mapping (^. 

/' = ^7/ (8) 

At linear order in 7, which will be valid for sufficiently weak shear, one can perform a first order Taylor 
expansion of the RHS of (0) and Sj becomes the differential operator 

5^ = 1 - ^aMaijTidj (9) 

where dj = d/drj. This operator is rather similar to a rotation operator. An important question is what 
is the domain of validity of (^. The answer depends on the content of the image to which it is applied. 
For an image containing information only at spatial frequencies below some upper limit fcmax this will be 
a good approximation provided r <C l/(7^max): so for finite shear (|^) will not apply for spatial frequencies 
^ l/(7r). The limit here is that combination of frequency and distance from the origin that a shear of 
strength 7 corresponds to a local translation of about an inverse wavenumber. We would however expect 
to correctly describe the effect of shear on the low frequency behavior of an image, in the sense that 
if applied to an image which has had the high frequency information removed, then the result will be 
essentially identical to the low-frequency content of an exactly sheared image. 



1.2. Second Moment Shear Estimators for Perfect Seeing 



A pure shear will cause a circular object to become elliptical, and will change the ellipticities of 
non-circular objects. Following Valdes, Jarvis, fc Tyson (1983| ) a natural choice of statistic to measure such 
a distortion is the second central moment or quadrupole moment 



qir. 



d'^r rirraf[r) 



(10) 



where we will assume that the total flux, which is unaffected by a pure shear, has been normalized such 
that / (Pr f{r) = 1, and where the origin of coordinates has been taken such that the dipole moment 
di = J (Pr rifir) vanishes. The triplet coefficients of the symmetric matrix qim are 



1 

-MA^Jq^J = 



{qxx + qyy)/'^ 
Qxy 



(11) 



The first component go is a measure of the size, or area, of the object, while the latter two are a 
measure of the eccentricity of the object — they vanish for a circular object — which we will refer to as the 
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'polarization'. We can compute how the quadrupolc moment is affected by a shear by applying (||) to / in 
(O) to find 



where we have integrated by parts and invoked the symmetry [M^ij = Maji) and tracelessness (M, 
of the matrices Mi, M2. With — MsimqE etc., and using(0), we find 

5qA = ^Malmiq'ljn ^ ffi™) = laM AlmMaliM Birnqs = 27/3((5^o9/3 + (^A/J^o) 



(12) 
= 0) 

(13) 
(14) 



or equivalently 

q'a=qa + 27a 90 

9o = 90 + 27Q(7a 

A weak gravitational shear therefore causes a change in the the polarization 5qa — 2qQja, which is 
proportional to the area qo, and it also induces a change in the area 6qo which is proportional to the 
eccentricity. For intrinsically randomly oriented galaxies the average intrinsic polarization (qa) vanishes by 
symmetry, so (q'^) — 2(go)7Q and (qq) = (qg) and therefore 



7a - {qDmq'o) 



(15) 



which which one can use to estimate the shear on a patch of the sky by replacing the averaging operator 
(. . .) by a summation: ■ ■ ■ f^- 



This type of shear estimator was first introduced and used by Valdes, Jarvis, & Tyson (1983) in their 
search for large scale shear. The averages of qa and qo here are heavily weighted towards larger galaxies, 
which is not ideal. More commonly, what has been done (e.g. Tyson et al. 1984 ) is to normalize the 
polarization by the trace of the second moments and define the 'ellipticity vector' 



Bq = qa/qo 

which depends only on the galaxy shape, and whose expectation value is 

' qa + 2jaqo' 



(e'a) = ( '"t^^"'° ) =^ 27. - 2^,{qaq^/q'o) = ^iM - (e^)/2) 

\ qo + 7/39/3 / 



(16) 



(17) 



and where we have kept only terms up to linear order in the shear ( Kaiser, Squires, fc Broadhurst 1995 
hereafter KSB). 



An alternative (Bonnet & Mellier 1995) is to normalize the second moments by the square-root of the 
determinant of qim rather than by the trace: 



qa 



\/\qirn\ \/ql~ qaqa 

which, again at first order in 7, has expectation value 



(er> = 27a^v^r7m^) 



(18) 



(19) 



Yet another possibility is to use only the information in the position angle 9 = (tan (72/91 )/2 (Kochanek 



1990| ), or equivalently in the unit shear vector |ea|, and there are numerous other similar statistics one could 



use. 
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These formulae for the response of the polarization statistics should be used with caution. While 
formally correct, the averages here must be taken in an unweighted manner over all galaxies. This is 
generally neither possible nor particularly desirable, as there are detection limits, and one would ideally 
like to estimate the shear as some optimally weighted combination of polarizations of galaxies of different 
fluxes and sizes, and these simple relations no longer hold. We will return to this later. First, however, let 
us consider the effect on these estimators of a finite point spread function, which is well illustrated by the 
simple estimator (|l5|). 



1.3. Second Moment Shear Estimators for Finite Seeing 

In real observations, some or all of atmospheric turbulence, optical aberrations; aperture size; guiding 
or registration errors; atmospheric dispersion; finite pixel size; scattering etc. will combine to give an 
observed surface brightness 

/o(r) = / d\' g{r')f{r -r') = g®f (20) 



where g{r) is the point spread function (PSF). The combined effect of gravitational shearing followed by 
instrumental and atmospheric seeing is the transformation 

fo^g^s.j (21) 

Circular seeing will tend to reduce the ellipticity while departures from circularity will introduce an 
artificial polarization. To make accurate shear measurements we need to correct for the latter and calibrate 
the former. For estimators like (p^, which are computed from moments qim as defined in (p^, the effect of 



the PSF is rather simple since, as noted by Valdes, Jarvis, & Tyson (1983), the second central moment of a 



convolution of two normalized functions is just the sum of the second central moments: 

qirnifo) = qimif) + Qlniig) (22) 

and this additive property is shared by the independent components qa- Thus one can recover the second 
moments (?a(/) that would be measured by a large perfect telescope in space from the observed moments 
QaHo) simply by subtracting the moments of the PSF qA{.f) = <lA{fo) — lA{g)- These may be measured 
from shapes of foreground stars in the image, or, in the case of diffraction limited seeing, computed from 
the telescope design. In terms of observed moments, the shear estimator (|l5| ) becomes 



where superscript * denotes the value for a stellar object. This procedure — using (|2^) to correct the 
measured moments of the PSF — simply and rather elegantly compensates for both the circularizing and 
distorting effects of realistic seeing. 



1.4. Weighted Moment Shear Estimators 

Unfortunately, while very simple to analyze, the shear estimator constructed from the second moments 
as defined in is not at all useful when applied to real data. For one thing, photon counting noise 
introduces an uncertainty in the moments which diverges as the square of the radius out to which one 
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integrates, and the effect of neigliboring objects will similarly grossly corrupt the signal. There is also the 
problem that for the kinds of PSFs that arise in real telescopes, the second moment is not well defined. To 
obtain a practical estimator, it is necessary to truncate the integral in dlQ). This can be done in various 



ways; the approach implemented in the FOCAS software package ( Jarvis fc Tyson 1981 ) and also the 



Sextractor package ( Berlin fc Arnouts 1996[ ) is to truncate the integral at some isophotal threshold /dit and 



compute moments of the non-linear function F(fo) = 0(/crit — /o)/oj where 6(/) is the Heaviside function. 
In fact, isophotal moments are most commonly computed from an image which has been smoothed with 
some kernel — usually an approximate model of the PSF itself — in order that the isophotal boundary be 
well defined. An alternative, as advocated by Bonnet fc Mellier (1995| ) and by KSB, is to limit the range of 



integration with a user supplied weight function w(r) and define 

qim oc / (fr w{r)rir.mfo{r). (24) 



Another possibility is to define a polarization vector in terms of the second derivatives of a smoothed image 

fs = w® J as 

qim « didm / fs (25) 



evaluated at the peak of fs- For a Gaussian weight function w{r), however, this is essentially equivalent 
to the weighted second moment (|2^). As we shall see, the weighted moment statistics, being linear in the 
surface brightness, offer significant advantages over the isophotal threshold method, but in either case, the 
simple relation (^2|) between observed and intrinsic second moments no longer holds, and compensating for 
the effects of the PSF becomes considerably more complicated than equation (p3|). 

A partial solution to this problem has been offered by KSB, who computed the response of shear 
estimators like ( p^ to an anisotropy of the PSF under that assumption that this can be modeled as the 
convolution of a circular PSF .gciic('') with some compact but possibly highly anisotropic function k{r). 
This would be a reasonable approximation for example for the case of atmospheric seeing in the presence 
of small amplitude guiding errors. They found that such a PSF anisotropy would introduce an artificial 
elli pticity Scg — P^g{fo)p f3 where pfj = Mpim J (Pr rirmk{r) is the unweighted polarization of k{r) (though 
see Hoekstra et al. (1998 ) who corrected a minor error in the analysis). The 'smear polarizability' P™ is a 



combination of weighted moments of /o, and is essentially a measure of the inverse area of the object. An 
interesting feature of this analysis is that only the second moment of the convolving kernel appears here; all 
other details of k{r) are irrelevant, and it is relatively straightforward to determine Pa from observed stars, 
set up a model for how this two- vector field Pa{r) varies across the field, and then correct, at linear order 
at least, the ellipticities to what would have been measured by a telescope with a perfectly circular PSF. 

KSB also computed the response of the ellipticity to a shear applied to /o after smearing with the PSF 
and found Sca — P2pifo)li3iQ) where Pj^ is another combination of moments of /□ . This of course is not 
what one wants, as one really needs the response to a shear applied before smearing with the PSF, and they 
suggested that one use deep images from the Hubble Space Telescope to empirically deduce a correction for 



finite seeing. A related approach, suggested by Wilson, Cole, fc Frenk (1996 ) is to iteratively deconvolve 



the images, apply a shear and then re-convolve and again use the change in the polarization with applied 
shear to calibrate the relation between 7 and the polarization measured from the original images. 

Luppino & Kaiser (1997) have used a somewhat different approach. They noted that the real operation 
( pl| ) can be written as 5 18) S^f = S-y{{S~^g) (8) /) i.e. applying a shear before smearing is equivalent to 
smearing with an anti-sheared PSF S~^g and then shearing. Now if the PSF is Gaussian, applying a weak 
shear to it is precisely equivalent to smearing it with another small but anisotropic Gaussian, so the effect of 
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this can therefore be computed using the smear polarizabihty of KSB, and it then follows that for a nearly 
circular Gaussian PSF, the operation ( [2l| ) will cause a response 

ea ^ = e„ + Se^ = e„ + P2^70 + P^]}pp{S-^g) (26) 

where pp{S^^g) is the unweighted polarization of the anti-sheared PSF and is of first order in 7. One 
can most easily infer the value of pi3{S^^g) from the values of the shear and smear polarizabilities for a 
stellar object: shearing a point source has no effect, so for a star we must have Pa = —{P'^{g)IP^™'{g))la 
where we have suppressed the indices on the stellar polarizabilities — for a nearly circular PSF these are 
approximately diagonal — and hence the net effect of a real shear in this approximation is 

•Jeo = {PlpUo) - {P'{9)lP'''\9))P^0{U))lp. (27) 

This is nice, as it expressed the linear response of the polarization Bq, to a shear entirely in terms of the 
observables /o and but it rests somewhat shakily on the assumption that shearing a realistic PSF can 
indeed be modeled as smearing it with some compact kernel. For a Gaussian PSF this is exact, but this is 
a rather special, and unfortunately unrealistic, case. 

One indication that (p7| ) cannot apply for a general PSF comes from considering the factor 
{P'^ig)/ P^^{g))- At no point have we specified the scale of the weighting function wir), so this factor must 
be invariant to choice of scale length. For a Gaussian PSF this is indeed the case, but for a PSF generated 
by atmospheric turbulence for instance this is not the case, and ( p7| ) is then inconsistent. Similarly, the 
factor qaid) / P^^{g) appearing in the KSB correction for PSF anisotropy is, in general, dependent on the 



scale of the window used to measure the PSF properties. Hoekstra et al. (1998) have found from analysis 



of simulated images that for the HST WFPC2 instrument one can adjust the scale of the weight function 



for the stars to render the calibration (27) and the KSB anisotropy correction reasonably accurate, but it is 
not clear that this will apply in general. Indeed, for diffraction limited seeing, the inadequacy of the KSB 
formalism has a deeper root. While the assumption that the real PSF can be modeled as a convolution of 
a circular PSF with some compact kernel k{r) may be a good approximation for atmospheric turbulence 
seeing with small amplitude guiding errors and such-like, as we shall see, this is not the case in general. 

The current situation is therefore somewhat unsatisfactory. In this paper we will develop an improved 
method of shear estimation which does not suffer from the inadequacies noted above and works for a PSF 
of essentially arbitrary form. The layout of the paper is as follows: In ^ we construct an operator which 
generalizes (^) to finite size PSF and which generates the effect of a gravitational shear on the observed 
(i.e. post-seeing) surface brightness Jo- We first show that, quite generally, the effect on fo = g <E) f oi 
a shear applied to / is equivalent to a shear applied directly to /□ plus a 'commutator' term which is a 
convolution of fo with a kernel jaHa{r) where Ha{r) may be computed from the PSF g. We explore the 
properties of this kernel for various types and PSF. The kernel Ha involves In 5 (where the tilde denotes 
the Fourier transform), which appears, particularly in the case of diffraction limited seeing, to be formally 
ill defined since g vanishes at finite radius. We show however, that the operator which generates the effect 
of a shear on a filtered image which has had frequencies close to the diffraction limit attenuated does not 
suffer from this problem. In §^ we specialize to weighted second moments, and compute their response to 



shear, both for individual objects {3.1 and for the population of galaxies of a given flux, size and shape § p.2| . 
In §3^ we show how to optimally combine estimates of the shear from galaxies of various different types. 
The method is tested using simulated mock images. In §^ we discuss the correction for PSF anisotropy, 
and propose a new technique involving measuring shapes from images which have been convolved with a 
re-circularizing PSF. We draw attention to a hitherto ignored noise related bias and show how this can be 



- 8 - 



analyzed and corrected for. In §|5| we summarize the main results, and outline how they can be applied to 
real data. In the analysis here we will draw heavily on the properties of real physical PSFs, and we include 
as an appendix a brief review of basic PSF theory and discussion of PSF properties as they relate to weak 
lensing observations. 



2. Finite Resolution Shear Operator 

We now consider how the shear operator is modified by finite resolution arising either in the 
atmosphere, telescope optics or detector. Let the unperturbed surface brightness (i.e. that which would 
have been observed in the absence of lensing) be 

/o = 3 ® / (28) 

so the perturbed surface brightness is 

/o = 5®^7/ (29) 
Fourier transforming, and using the result that, at linear order in 7, applying a shear in real space is 
equivalent to applying a shear of the opposite sign in Fourier space (i.e. if a = S-yb then d = S-^b, where 
F{k) = J (Pr F{r)e'^-'"), we have 

/^=g5_,/ = /o-g55^(/o/ff) (30) 

where 5S^ = 5*^ — 1. Since 5S^ is a 1st order differential operator we have 5Sy{fo/g) — g^^SS~ffo~g^^ fo^S^g, 
and g~^dSjg — SS-ylng. Consequently, to first order in 7 

f:,^ fo-SSjo + foSSylng (31) 

or, in real space, 

= /o + SSyfo - (SS^h) (g, /o (32) 

where h is the inverse transform of /i = In^, i.e. the logarithm of the optical transfer function (OTF). 
Invoking the definition of the shear operator (H), we have 

/o = /o - laMaij{ridjfo ~ (ridjh) (g) /o). (33) 

Note that the PSF is a real function so g{—k) — g*{k), and this symmetry is shared by Ing, so h is also real. 

The finite resolution shear operator ( |33| ) is then — S-yfo — {SS^h) (g) that is, it is the regular shear 
operator S-y applied to the post-seeing image /o plus a 'commutator term' g (g) S-yf — Sry{g iSi f) which is a 
correction for finite PSF size and which is a convolution of /o with a kernel ^aHa{r) — ^aMdjridjh = SSyh 
which one can compute from the PSF g. This seems quite promising; the effect of the first term on the 
polarization of an object is just 7^ times the KSB post-seeing shear polarizability. The response of the 
polarization to the second term should also be calculable; if the kernel is very compact then the response 
will be given by the KSB linearized smear polarizability, but even if it is not, ( ^ ) should still allow one to 
compute the polarization response since it expresses the change in /o, and therefore in the polarization, or 
indeed any other statistic computable from directly in terms of the observed surface brightness /o itself. 
The finite resolution shear operator involves the function h{r) which is the transform of the logarithm of 
the OTF. Since the OTF becomes exponentially small or may vanish, two questions immediately arise: Is 
the function h{r) mathematically well defined? and can it be reliably computed from PSFs measured from 
real stellar images? To address these questions we now explore the form of this 'commutator kernel' Ha (r) 
for various types of PSF. 
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2.1. Gaussian Ellipsoid PSF 

Consider first the simple though unreahstic case of a Gaussian ellipsoid PSF: g(r) = 
(27r)^^ |m|^^/^ exp(— rim~-^rj/2), where = is the matrix of central second moments. In 

this case, In^(fc) = —kimijkj /2, and so dS'yhig(k) = jaMaim'nT'mjkikj and on transforming this we find 
that the kernel is the operator Ha = Maim'mmjdidj (in which case the function Ha{r) can be realized, for 
example, as the contraction of the constant matrix jaMaimmmj with the limit as cr — > of the matrix of 
second partial derivatives of a Gaussian ball (27ro'^)^^ exp(— r^/2o'^)), and therefore 

/o = /o - laMaijiridjfo + mudidjfo)- (34) 

So, for a Gaussian PSF, the finite resolution shear operator is well defined and is a purely local differential 
operator. Gaussian PSF's are, however unphysical, and do not arise in real instruments. 



2.2. Atmospheric Turbulence 

Now consider atmospheric turbulence limited seeing. As reviewed in appendix ^ in that case 
g{k) — exp{—S{kDX/2TT)/2) where for fully developed Kolmogorov turbulence the structure function 
is S{r) = 6.88{r/rQ)^/^ with rg the Fried length, so in this case ft, cx Ing cx —k^^^. The commutator 
kernel therefore involves the transform of this power-law, which diverges strongly at high k, so how do we 
make sense of this? Dimensional analysis would suggest that h{r) = J d?k k^/'^e^'^'^ be a power law with 
h{r) (X r~^^l^ . The same argument, however, applied to a Gaussian would suggest h{f) oc r~^, which we 
know to be false as we have just shown that in that case h[r) is just the second derivative of a (5-function 
and has no extended tail. To clarify the situation, and to verify the validity of the power-law form for h(r) 
for atmospheric seeing, consider the function h(r\ R) which is the transform of k^^^ times an exponential 
cut-off function exp{—kR), that is 

ft(r; i?) cx y d2/cA;5/3g-fei?eife r (35) 

so h{r) is the limit as i? ^ of h(r] R). If we write this as an integral with respect to rescaled wavenumber 
y — kR it becomes clear that /i(r; R) obeys a self-similar scaling with respect to choice of R and can be 
expressed in terms of some universal function F{y) such that h{r]R) = R^^'^/'^F{r / R). If we postulate 
that h{r; R) tends to some i?-independent limit for finite r as i? — > then that limit must be a power law 
with F{z) cx z~^'^/'^ so h{r) must be proportional to r~^^/'^. This argument does not indicate the coefficient 
multiplying the power law which, for a Gaussian, happens to vanish. The value of h at the origin is 
h{Q-R) = 27ri?-"/3r(8/3), i.e. on the order of the value of the r ^ R power law asymptote extrapolated to 
r ^ R (note that for the Gaussian the analogous gamma function is not defined). A numerical integration 
for various values of R is shown in figure |l| This shows that as we decrease R, the function h{r; R) does 
indeed tend to a t-^h/^ power law with finite non-zero i?-independent amplitude, but that the power law 
breaks (with a change of sign) at r ^ i? and becomes asymptotically flat for r ^ R. For finite R this 
function is well characterized as a positive 'softened (5-function' core with width ~ R, central value ~ i?^^^/"^ 
and therefore with weight proportional to R~^/^, surrounded by a negative power law halo with h cx r~^^^^. 
In the limit R 0, the core shrinks and becomes negligible, leaving only the power law halo. This power 
law has the same slope as the well known large angle g{r) cx r~^^^^ form of the PSF, but the PSF departs 
from this law at angular scale Vg ~ A/rp corresponding to the Fried length, whereas h(r) is a perfect power 
law and shows no features at the Fried scale. 
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Fig. 1. — Numerical calculation of h{r: R) = J cPk k^/^ exp( 
by factors of 2. Also shown is a pure r~^^/^ power law. 



—kR + ik ■ r) for a range of R values increasing 
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The function h(r) diverges strongly at the origin and the same is true of the kernel Ha = Maijridjh, 
so Ha = — ll/3r~^^/'^{cos2(^, sin2(^}. This does not give rise to any inconsistency with ( |3^ ) however. The 
observed sky /o is coherent on the scale of the PSF rg, so in computing the contribution to the commutator 
term 

Sfo = laMa^, J <fr' r[ry )-^'/^Ur - r') (36) 

at small r' <C ^o, we can perform a Taylor expansion of and we find that the first non- vanishing term 
is the second order term 5fo ^ {didjfo)r J d^r' r^r^ (r')"^^/'^ which has no physical divergence at r' = 0. 
Somewhat unfortunately perhaps, the same line of argument shows that the commutator term cannot in 
this case be assumed to be a convolution with some compact function k{r). A necessary condition for this 
to be valid is that the unweighted second moment of the kernel Pa — Maimlp J ^ifmHp{r) should be 
well defined and tend to a finite limit within some small radius ^ r^, the characteristic width of the PSF. 
But this is the same integral as above, which does not tend to any well defined value, but rather diverges 
as the 1/3 power of the upper limit on the integration radius. Thus while the h{r) cx r~^^^^ asymptote 
is quite steep, it is not sufficiently steep to render the pa value well defined. This strictly invalidates the 



approximation of Luppino & Kaiser (1997) for turbulence limited seeing 



The Kolmogorov law is only expected to apply over a finite range of scales. The structure function will 
fall below the r^^^ form at the 'outer scale', which recent measurements at La Silla suggest to be typically 



on the order of 20m (Martin et al. (1998), though see also the review of earlier results in Avila et al 
(1997 )). Fast guiding (often referred to as 'tip-tilt' correction) would also eff'ectively reduce the structure 
function at these frequencies, and such effects act in a similar manner to the simple exponential cut-off' we 
have assumed. To compute the OFT properly, one should include the effect of the aperture. The detailed 
form of h{r) at very small r will be sensitive to the detailed form of the outer scale cut-off and/or aperture, 
but provided these lie at scales much greater than the Fried length (which is the case for large aperture 
telescopes at good sites) the effect on the commutator term should be almost independent of the cut-off 
because any signal at the relevant spatial frequencies will have been attenuated by a large factor. There 
will also be deviations from Kolmogorov structure function law at small separation; diffusion will damp out 
fine scale turbulence, and mirror roughness will add additional high spatial frequency phase errors. These 
effects will modify the h oc r~^^/^ halo at large r, just as they modify the large-angle g oc r~^^/^ form of the 
PSF. These effects may profoundly inffuence the behavior of h, Ha at large angle, but have little impact 
on the type of shape statistics considered here, where the large-angle contribution to the polarization is 
suppressed by the weight function or the isophotal cut-off. 



2.3. DifTraction Limited Seeing 

Let us now consider diffraction limited observations. In this case the OTF g falls to zero at finite 
spatial frequency — the diffraction limit — and so In g diverges as one approaches the diffraction limit 
from below and is not defined for higher frequencies. It is easy to understand how this divergence arises 
physically. As noted above, applying a weak shear in real space is equivalent to applying a shear of the 
opposite sign in Fourier space. Thus the information in some Fourier mode of the sheared image comes 
from a slightly displaced mode in the unsheared image. If the OTF g is finite and continuous, as is the 
case for atmospheric turbulence limited seeing, then the image /o contains all of the information required 
to predict /q. For diffraction limited seeing however, and for observations through a narrow band filter, 
the OTF has a well defined edge, so the information contained in /□ for spatial frequencies just inside the 
cut-off may lie outside the cut-off in and so will be missing. For a finite band pass the form of the OTF 
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will be modified, but must still fall to zero at the diffraction limit for the highest frequencies passed by the 
filter and In 5 is still formally divergent. 

Thus, in general, from knowledge of /o = g ^ /, it is strictly speaking not possible to say how /o 
would change in response to a small but finite shear applied before seeing; there are Fourier modes within a 
distance Sk 7fcmax of the diffraction limit that one cannot predict. As an extreme example, imagine we 
have a pure sinusoidal ripple on the sky which lies just outside the diffraction limit and which is therefore 
invisible. Applying an appropriate shear can bring that mode inside the limit and the ripple will appear 
as if from nowhere. As applied to real signals, however, this formally divergent behavior does not present 
a serious problem. First, if we observe galaxies of overall extent re, which is typically not much larger 
than the seeing disk, then the transform must vary smoothly with coherence scale Sk ^ l/rc; which, for 
sufficiently small 7 will greatly exceed 7A;inax; this rules out the possibility of isolated spikes lurking just 
beyond the diffraction limit as in the example. Also, for filled aperture optical telescopes the marginal 
modes axe quite strongly attenuated, and, for any finite measurement error from e.g. photon counting 
statistics, will contain very little information. The situation here is similar to that encountered in analyzing 
atmospheric turbulence PSF; there, while the very small scale details of h{r) are sensitive to the aperture or 
outer scale cut-off, when applied to real data they have essentially no effect. Similarly, if one can generate a 
function which accurately coincides with h = log g where g{k) exceeds some small value, but which tapers 
off smoothly at larger frequencies (rather than diverging at fcmax), then this should have a well defined 
transform and should give what is, for all practical purposes, a good approximation to the true finite 
resolution shear operator. 

One way to explicitly remove the divergence is to compute shapes from an image which has had the 
marginally detectable modes attenuated. If one rc-convolvcs the observed field /o with some filter function 
g^r) to make a smoothed image /s = 5^ ® /o then in Fourier space we have 



so provided g^ falls off at least as fast as g as one approaches the diffraction limit this operator is well 
defined. In real space the corresponding operator is 



so now the commutator term is the convolution of /o with Ha = Maijg^ (S) {ridjh). In principle one can 
design g^ such that fs is essentially identical to fo', simply set g^ to unity for all modes within some tiny 
distance of the diffraction limit and zero otherwise. However, this may not be a very good idea; the sharp 
edge of such a filter function in fc-space will result in ringing in real space, both in the filter .9^(r) and 
especially in Ha{r). These extended wings have no effect on real signal, but will couple to the incoherent 
noise in the images, whose spectrum does not fall to zero as one approaches the diffraction limit. To 
ameliorate these problems one can choose g"^ to have a soft roll-off as one approaches fcmax to counteract the 
divergence of 1/g. One simple option is to set g^ = g; i.e. to re-convolve with the PSF itself. In this case 
Ha{r) will be about as compact as the PSF, and we then have 



fl = h- g^SSjo + fog^g-'SS^g 



(37) 



/s = fs - laMaijig^ ® {ridjfo) - g^ ® (ndjh) ® /o) 



(38) 



fs 



fs - laMaijig ® {Tidjfo) - {fo ® {ridjg)) 



(39) 




(40) 
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To summarize, we have shown in ( p^ ) that the effect on the observed sky /□ of a weak shear appUed 
before seeing can be written as a shear apphed after seeing plus a convolution with some kernel which 
is the sheared transform of the logarithm of the OTF, which can be computed from the PSF. We have 
explored the form of this kernel for both Gaussian and more realistic models for the PSF. For atmospheric 
turbulence limited seeing the kernel is a power law H ~ j.-ii/3 diffraction limited seeing the shear 
operator appears formally ill-defined but this is not a serious problem when applied to real data and that 
one can, for example, compute the operator for a slightly smoothed sky fs=g'^® fo in a divergence free 
manner. We have presented the shear operator for the case where g'^ — g. This could potentially be used to 
compute the response of the shape statistics measured by the FOCAS and/or Sextractor packages, as these 
are measured from just such a re-convolved image, but we will not pursue this here. 



3. Weighted Moment Shear Estimators 

We now specialize to weighted quadrupolc moments as defined in (|2^ ) and compute how these respond 
to shear. We first compute the response of the moments of an individual object, and we then compute 
the conditional mean response for a population of objects having given flux, size etc. This will allow us to 
compute an optimal weight function for combining shear estimates from galaxies of a different fluxes, sizes 
and shapes. 



3.1. Response of Weighted Moments for Individual Objects 

Consider again for illustration the case of a Gaussian ellipsoid PSF g{r) cx exp{—rimJ'j^rj/2) and 
moments qim = / d'^r foir)w{r)rirm. From (|3|), ^) we have ql^ = qim + Sqi„i with 



Sqim = -JaMatj J (fr wrir,n{ridj fo + mipdpdjfo) (41) 

Integrating by parts to replace derivatives of fo with derivatives of wriVm we find the linear response of 
qji = ^Mjiimqim to a shear can be written 

SqA = Pap1(^ (42) 

with 'shear polarizability' 

PAf3^ I d\VAp{r)fo[r) (43) 



and where ^ 

'PAp{r) = -MAimMi3ij[ridjWrirrn - rriipdpdjwrirm] (44) 

For the special case of w = 1, i.e. unweighted moments, wc find Vap — 5ap{riri — ma) and hence 
5qa — ^aiqu — fn-ii) — 2ja{qQ ~ ™o) in accord with the result obtained in the Introduction. In this case, 
Paf3 is a combination of zeroth and second moments of fo- 

For a general PSF and for moments measured from a filtered field fs — g'' fo as a weighted moment 

qA^\MA^^jd-r^ir)nrMr) (45) 
we find using (^3|) that the response 5qA can be cast in the same form, but now with 

'PAp{r) = ]^MAimMp,j[ri{g^ © {djLunr^)) - {nh) ® g^ ® (a^wr/r™)] (46) 
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where we have defined the correlation operator such that (a © 6)^ = / d?r' a{r')h{r' + r). If the moments 
are measured directly from the unfiltered images /o then one can replace g'^{r) with a Dirac (5-function to 
obtain ^ 

VApir) = -MAimMp^jlnd-jWrirm + (hu) © (djwnrm)] (47) 

whereas for the special case — g 

VApir) = ^MAimM0,j[2ri{g © dj{wrir„,)) -g® {ridj{wrir„,))]. (48) 

The function Vapir) is shown in figure ^jfor a turbulence limited PSF and for a Gaussian window function 
wir). Note that ( ^6|) and (^8|) are well defined continuous functions even in the limit that the weight 
function becomes arbitrarily small; i.e. w{r) — > 5{r). 




Fig. 2. — The panels on the left show the pair of functions wi, W2 with 'Wa{r) — \Maimw{r)rirra which 
when multiplied by /s and integrated give the polarization statistic qa- The four panels on the right show 
the components of the polarizability kernel Va^ir) which when multiplied by /o and integrated yields the 
polarizability Pa/B- The PSF was computed from a turbulence limited OTF g = exp(— 0.5(fcr*)^/'^), and the 
smoothing kernel was w{r) = exp(— 0.5(r/r*)^) with scale length r* equal to one eighth of the box side. 

Equation (^2|), along with the appropriate expression for PA^(r) tells us how the polarization statistic 
for an individual object formed from weighted quadrupole moments responds to a gravitational shear. This 
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is an essential ingredient in calibrating the shear-polarization relation for a population of galaxies. As we 
shall see in the next section, there are some subtleties involved, but for now we note that if we simply 
average (|2|) over all galaxies on a patch of sky we have 

{q'a) - ila) - {Pcflhp (49) 

SO an estimate of the net shear is given by 

%^{Pc.p)-'{{q'0)-{qp)) (50) 



Which is the generalization of (^3|) to weighted moments. The term (q'p) in ( pO[ ) is the averaged observed 
polarization. The term (qa) is the mean polarization generated by anisotropy of the PSF, and we will show 
how this can be dealt with below. 

Ideally, in averaging polarizations, one should apply weight proportional to the square of the signal to 
noise ratio, which one would expect to be a function of the flux, size, eccentricity etc. of the objects. If 
the shape is measured with some fairly compact window function w{r), then the total flux, which may be 
dominated by the profile of the object at considerably larger radius, is probably not ideal and one will likely 
obtain better performance if one takes the weight function to be a function of go, = qaqa and a weighted 
flux 

F= j (frw{r)fs{r). (51) 



The response of F can be computed in much the same way as qA, and we find from (jS^) 5F — Ra^a with 
i?a = j (fr Uafoir) with 



or, for the case = g 



n^{r) = M„y[r,(5-t © d^w) - {r,h) © 5^.^^)] (52) 



7^„(r) = Maij[2ri{g ® d^w) -g® djwn)]. (53) 



In the foregoing we have implicitly assumed that applying a shear does not affect the location of 
an object. This is not necessarily the case. If objects are detected as peaks of the surface brightness /□ 
smoothed with some detection filter Wd, that is as peaks of fd — Wd ® foi then for an object which in 
the absence of shear lies at the origin, we have di{{)) = di{wd ® /o)o — J d^i" diWd{r)fo{r) = while after 
applying a shear we have di{0) = ^aM^aim J (Pr diWd[ridmfo — {ridmh) (8) /o]- This will not, in general, 
vanish, implying that the peak location will have shifted, and consequently the central second moments 
should be measured about the shifted peak location, whereas in the above formulae we have computed the 
change in the moments without taking the shift into account. One could incorporate this effect, but at the 
expense of considerable complication of the results. There is some reason to think that this effect is rather 
weak. In particular, if the galaxy is symmetric under rotation by 180 degrees, so fo{f) = /o(— then the 
shift in the centroid vanishes. In general this is not the case, and the formulae above should be considered 
only an approximation. 



3.2. Response of the Population 

Equation ( pO| ) above gives a properly calibrated estimate of the gravitational shear. It is, however, less 
than ideal as the polarization average is taken over all galaxies with equal weight. This is neither desirable 
nor is it achievable in practice due to selection limits, and what one would rather have is an expression 
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for the average induced polarization for all galaxies in some cell of flux, size and shape space, which we 
parameterize by F, qo, and = qaqa- One can then average appropriately weighted combinations of the 
average shear estimate for each cell. The mean induced polarization for such a cell depends not only on the 
polarizabilities of the objects contained therein, but also on the gradient of the mean density of objects as 
a function of the photometric parameters F, qo, qa- Consider a slice through this 4-space at constant F, qo. 
A shear will induce a general flow of particles in this space in the direction qa — Ja- The mean polarisation 
for a cell in F — qo — q^ is the average around an annulus in q^ space, and depends quite sensitively on 
the local slope of the distribution of particles in \qa. These factors can have a profound influence on the 
weighting scheme; for a distribution which is flat near the origin in qa space, like a Gaussian for example, 
nearly circular objects have no response and should therefore receive no weight. For a randomly oriented 
distribution of circular disk galaxies, in contrast, the distribution in Qa-space has a cusp at the origin; the 
response becomes asymptotically infinite, and these objects dominate the optimally weighted combination. 
These examples are both idealized, but underline the importance of computing the population response in 
order to obtain optimal signal to noise. 

Let us first compute the conditional mean polarization for galaxies of a given flux and size: (qa) F,qo ■ 
The mapping of the photometric parameters F, go, Qa is 

F' = F + Rp-ip 

q'o^qo + Popip (54) 

so we need to consider the distribution of galaxies in F, qo, qa, Ra, Pna, Pap with lensed and unlcnsed 
distribution functions related by 

n'{F', q'o, q'^, R'^, P^,^, P'^p)dF' dq'^cp q' CP R(f P^^d^ P' - n{F, qo, qa, Ra, Poa, Pap)dFdqod^qd^ Rd^ Pod^ P. (55) 

Multiplying by W{F' , q'o)q'a, where W is some arbitrary function, and integrating over all variables we have 

J dF'dq'od^q'd^Rd^P[^d'^P'n'W{F',q'Q)q'^= J dFdqod^qd^ Rd^Pod'^PnW{F + 5F,qo + Sqo){qa + Sqa). (56) 

Now to zeroth order in 7 this vanishes because of statistical anisotropy of the unlensed population, so using 
(p4) for 6F etc. and performing a Taylor expansion of W{F + 5F, qo + Sqo) and integrating by parts we have 



JdF'... d^P' n'W{F', q'o)q'a ^ J dF . . . d-'P W{qo) (^nPap - ^o^9a^ - Rpqa 



dn 
'dF 



(57) 



To first order in 7 we can replace unprimed by primed quantities throughout the integral on the RHS since 
we need to compute this only to zeroth order accuracy. We now have a relation between the mean value of 
the observed polarization on the LHS and some other observable, again integrated over the distribution of 
observed galaxy properties (rather than of the unlensed parent distribution). Since W{qo) is arbitrary, and 
dropping primes, this implies 

J d\d^Rd^Pod^Pnqa =^pJ---J d\d^Rd^Pod^P (^nPap - ^opQa^^ - Rma^^ (58) 
or equivalently, that the conditional average polarization is 

{qa)F,qo = 7/3 



, I dn{Popqa) I dniRpqa) 

\Pap) - -■ 



n dqo n dF 



(59) 
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where n = n{F,qQ). Thus, as expected, the shear induced shift in the mean polarization for galaxies of a 
given size and flux differs from the mean of the shift 6qa = Paplp for an individual object. This is because 
a shear changes the weighted flux and size of an object in a way which is correlated with its shape. When 
we average the shear for galaxies in some cell in flux-size space we are averaging over galaxies which have 
been scattered in size and flux and we obtain a bias in the net polarization which depends on the gradients 
of the distribution function. 

We can generalize this to compute the response for galaxies of given flux size qo, and rotationally 
invariant shape parameter q^ = qaqa- To do this, we set q^ = qq^ with the unit polarization vector 
qa = {cosip, sin 1^9}, so (Pq = qdqdip. We then have, now for some arbitrary function W{F, qo, q^), 

j dF'dq'Qdq^'dip'(fR(fPod'^P'nW{F',q'Q,q^')q^ = J dFdqodq^dipd^Rd^Pod'^PnW{F+6F,qo+Sqo,q^+6q^){qa+Sqc) 

(60) 

where now n — n{F, qo,q'^, (p, Ra, Poa, Paf})- Using 5F etc. from (|5^ ) and 5q^ — 2qjjPjfi3^i3 we have 
or equivalently {qa)F.qo.q^ = Paplf3 with effective polarizability 



-D /D \ "2 dn{q,,P,jf3qa) I dn{Papqa) I dn{Rf3qa) , , 

Paf3 = {Pap) ^ ^ 7^ (d2) 

n oq'^ n oqa n ot 



where now n = n{F, qo,q ), and all averages are at fixed F, qo, q . 

The photometric parameters qa,qa appearing here are unnormalized. This is convenient for computing 
the linear response functions. It is, however, somewhat awkward here since the distribution function 
n{F, goj 9^) is highly skewed since go and q^ correlate very strongly with the flux. In computing the effective 
polarizability it is more convenient to work with rescaled variables q^ = qo/F, and q^' = q'^/F'^. The 
distribution function in rescaled variables is 

n'{F,q'o,q^')^FMF,Qo,q^)- (63) 

If we also re-scale the polarizabflities R'^ = Ra/F, P'j^p = Pap/F, p'j^^ — Pa^/F, re-express (^) entirely 
in terms of primed quantities and then drop the primes we flnd 



p _ ,p X 2 dnjqrjP^pqg) 1 dn(Popqa) ^ 1 
n oq'^ n oqo n 



^ d d , „ d 



dlnF dlnqo din q^ 



n{Rpqa) (64) 



where we have used the result that for any function X{F, q'^, q^') the partial derivative WRT F at constant 
qo, q^ is idX{F,q'o,q^')/dF)q„^g2 = dX/dF - {q'g/F)dX/dq'o - {2q^'/F)dX/dq^' to compute the term 
involving R^. The rather cumbersome expression (^) calibrates the relation between the shear and the 
mean polarization for galaxies in a small cell in flux-size-shape space. To compute it we need to bin galaxies 
in this space to obtain the mean density n and the various averages appearing here, and then perform 
the indicated partial differentiation. The form ( |6^ is somewhat inconvenient as the density n{F, qo, q^) is 
asymptotically constant as and one has to properly deal with the discontinuous derivative at this 

boundary. A computationally more convenient approach is to make one final transformation from q^ — s- g; 
since n{F, qo,q) — 2qn{F, qo, q^) falls to zero as q — > 0, and there is then no need for any special treatment 
of the derivatives at the boundary. With this transformation we have 



-D /D \ I dn(q^Pr,(3qa) 1 dn(Popqa) , q 

Pap = {Pap) ^ ^ \- - 

n oq n oqo n 



Odd 

1 - 



9 In F d In qo dlnq 



n{Rpqa) (65) 
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where now n = n(F, qQ,q). 

What about measurement noise? Let us assume that one has been given an image containing signal 
and measurement noise, and that one has detected objects, and measured quantities hke F, qj^. How would 
these photometric parameters change under the influence of a gravitational shear? The answer is given by 
(U), but with the understanding that Rp, Fo/3 etc. be the response functions one would measure in the 
absence of noise. This means that ( |65|) is also applicable with the same proviso. The major terms in ( |65| ) 
are however invariant of additive noise. The exceptions are the terms invoving (Popqa) and {Rpqa) which 
are quadratic in the sky surface brightness. The measured expectation values (Popqa) etc. therefore exceed 
the true values, but by an amount one can calculate from the known properties of the measurement noise. 
Another implicit assumption in the above analysis is that the objects are actually detected, which restricts 
applicability to objects which are detected at a reasonable level of significance. Aside from this, the results 
above should be applicable in the presence of measurement noise. To test these claims we have made 
extensive simulations with mock data, the details of which are described in appendix Figure || shows the 
results of one of these. The actual polarization agrees quite closely with the effective polarizability, even for 
very faint objects. While the differences between the effective polarizability and that for individual objects 
is not very large - typically on the order of 20% or so - the effective polarizability clearly describes the true 
response more faithfully. We shall now use (^5|) to construct a minimum variance weighting scheme for 
combining shear estimates. 



3.3. Optimal Weighting with Flux, Size and Shape 

Armed with the conditional mean polarization {qa)F.qo.q^ '^e can now compute an optimal weight (as 
a function of, for example, flux F, size qoi and eccentricity q^) for combining the estimates of the shear 
from galaxies of different types. Let us assume that one has measured fluxes etc. for a very large number 
of galaxies — from an entire survey say — and that from these data one has determined the mean number 
density of galaxies n{F, go, 9^) and also the various conditional averages and gradients that appear in (|6^). 
Now consider a relatively small spatial subsample of these galaxies and bin these into cells in F, qg , q^ space, 
and for each bin compute the occupation number N and the summed polarization ^ ■ 

A shear estimator for a cell which happens to have occupation number N is 

To simplify matters, let us neglect for now any anisotropy of the point spread function, in which case we 
can write Pap — PSajS where P = Pijr,/2 and we then have 

The expectation value for the variance in 7 (for a cell which happens to contain N galaxies) is 

(7') = (E l')/iN'P') = q'/NP' (68) 

where we have used {{J2la){J2la)) — since, in the hmit of weak shear, the galaxy polarizations 

are uncorrelated. As different cells give shear estimates whose fluctuations are mutually uncorrelated, the 
optimal way to combine the shear estimates from all the cells is to average them with weight per cell 
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Fig. 3. — The top panel shows the density of galaxies detected in the simulations described in appendix^ as 
a function of F, qq, \q\. Each of the sub-plots shows the density at fixed F with abscissa j^l and ordinate qq. 
These sub plots are arranged with flux increasing from left to right with flux F increasing by a factor 1.43 
at each step. All particles detected with significance level greater than 4-sigma are shown. The panel below 
this shows the density of particles weighted by the individual object polarizability, and the panel below that 
shows the density weighted by the effective polarizability given by (jSSj). The lower plot shows the density 
of objects weighted by the actual polarization (divided by the shear applied in the simulation, or 7 = 0.1 in 
this case). 
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2 

Wcoii oc 1/(7^) = NP jq^ to obtain a final optimized total shear estimate 

-total _ cells _ galaxies _ galaxies , , 

cells galaxies galaxies 

where Q = P/q. Thus the optimized cell weighting scheme corresponds to averaging the shear estimates 
from individual galaxies ^g^'^'^y = qa/P (for galaxies of given F, go 7 9^) with weights Wgaiaxy ^ P /q'- 

The variance in the total shear estimator is 

(7^) = ^^ = (E^?V^ (70) 

The quantity ^ is extensive with the number of galaxies and its value, per unit solid angle of sky, 
provides a useful figure of merit for weak lensing data. From a 2.75hr I-band integrations of solid angle 
dn = 0.165 square degrees at CFHT taken in good seeing (0".60 FWHM) we obtained J2Q'' - 4-7 x 10"^ 



( Kaiser et al. 1998 ) or ^Q^/<Kl ~ 2.85 x lO'^/sq degree, so with data of this quality, the statistical 
uncertainty in the net shear (per component) measured over one square degree would be around 
cr^ ~ (2 X 2.85 X 105)-i/2 = 1.32 x lO'^. 

This figure of merit allows one to tune the parameters of one's shape measurement scheme, such as the 
weight function scale size, in an unbiased and objective manner. The weighting scheme derived above is 
appropriate if the shear is independent of the measured flux etc. This is the case for lensing by low redshift 
clusters where for all relevant values of source redshift the sources are effectively 'at infinity' and the shear 
has saturated at its value for an infinitely distant source 7^. For high redshift lenses the shear will vary 
with source redshift, and if one has some, perhaps probabilistic, distance information at one's disposal, 
then the weighting scheme should be modified. Let us assume that the measured photometric properties pi 
of some object indicate it has a probability distribution to be at distance z oi p[z\pi); with high resolution 
spectroscopy this would be a delta- function at the measured redshift whereas with only broadband colors 
the conditional probability would be smeared out and perhaps multimodal. The conditional mean shear for 
this object, for a given a foreground lens, is proportional to the mean inverse critical critical surface density, 
and the optimal estimate for the shear at infinity is 

(l/Scrit(2))Qga 

- 00 ^ 1 galaxies , , 

Ecrit(oo) E (l/Scrit(^))2Q2 ^ ^ 

galaxies 



4. Correction for PSF Anisotropy 

In the foregoing we have computed how the shape polarization qa responds to a shear. This allows us 
to correctly calibrate the circularizing effect of seeing. In general, asymmetry of the PSF will also introduce 
spurious systematic shape polarization in (pO[), and it is crucial that this be measured and corrected 
for. As discussed in sQ, for unweighted second moments, the effect of PSF anisotropy is rather simple since 
the (flux normalized) final second moment qim is the sum of the intrinsic second moment and that of the 
instrument, so the anisotropic parts of the instrumental second moment qa — Mahnlim can be measured 
from stars, and subtracted from those observed. For weighted or isophotal moments things are more 
complicated, and depend on how the anisotropy is generated. KSB considered a simple model in which 
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the final PSF is the convolution of some circularly symmetric PSF with a small, but highly asymmetric 



anisotropizing kernel' k{r). In §4.1 we will further explore this 'convolution model', the attempts that 



have been made to implement this perturbative approach, and estimate the error in this method. In §4.2 



we show that the convolution model is quite inappropriate when applied to diffraction limited seeing and 



we develop a more general technique for correcting PSF anisotropy. Finally, in §4.3 we draw attention to 
a type of noise related bias that affects shear estimators, but which has hitherto been overlooked, and we 
show how this can be dealt with. 

The artificial polarization produced by PSF anisotropy is an effect which is present in the absence 
of any gravitational shear. Thus at leading order in 7 we can set 7 = 0. This effectively decouples 
the computation and correction of the PSF anisotropy from the shear-polarization calibration problem 
considered above. It also means that in this section we can assume that the intrinsic shapes of galaxies are 
statistically isotropic. 



4.1. The Convolution Model 

KSB computed the effect of a PSF anisotropy under the assumption that this can be modeled as the 
convolution of a perfect circular PSF with some compact kernel fc(r). This would be a good model, for 
instance, for observations primarily limited by atmospheric turbulence with seeing disk of radius Tg, but 
with very small guiding or registrations errors or optical aberrations with extent Sr ^ Vg, and the KSB 
analysis gives a correction to lowest order in Sr. In this model, the effect of 'switching on' the anisotropy is 
the transformation on the observed image, which is smooth on scale ~ r^, 

fo{r) fL{r) = J d\' k{r'),Ur - r') = /o(r) - kA.Ur) + h^^d.d.Ur) ~ h,jAd,dkIo{r) + . . . (72) 

where we have Taylor expanded fo{r — ?■') and where ki = j (Prrik{r), kij = J d?rrirjk(r) etc. Taking the 
kernel to be centered, we have h — 0, and to lowest non- vanishing order the effect on the polarization is 



d^r w^{r)f^{r) = + -k,, / d\ Ur)d,d,w^{r) (73) 



where Wa = Maimw{f)^ifm and we have integrated by parts. The induced stellar polarization is linear in 
kij and therefore scales as the square of the extent of the convolving kernel since kij ~ 5r'^ . Performing the 
decomposition kij = kAMAij we find that an individual galaxy's polarization qa will depend on both the 
trace ko and the trace- free parts ka of kij. The average induced polarization, however, only depends on ka 

and we have {qa) — kpP^ with 



PTp = T^Mai^Mp,, / (fr {fo)d^dj(w{r)rirra). (74) 



To lowest order in the PSF anisotropy we can use the KSB polarizability (Q) to infer k^ from the shapes 
of stars, and then correct the weighted second moments or the galaxies. At the same level of precision 



one can convolve one's image with a small kernel designed to nullify the anisotropy. Fischer & Tyson 



(1997 ) have presented a 3 x 3 smoothing kernel which does this. In the typical situation, the shapes are 
measured from an average of numerous images taken with a pattern of shifts, and which are co-registered 
to fractional pixel precision, and a simpler but equally effective approach is then to average over pairs of 
images which have been deliberately displaced from the true solution by the small displacement vector 
±Sr = y|fcI|{cos(6'),sin(6l)} with 9 ^ tan-i(fc2/fci)/2 + Tr/2. 
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What is the error in this hnearized approximation? To answer this we must consider higher terms in 
the expansion ([72|). The next order correction to qa involves kijk which, unhke the centroid ki cannot be 
set to zero, so for a given object the fractional error in the KSB correction is on the order of 6r/rQ ~ ^/e^- 
Galaxies, however, are randomly oriented on the sky, so the average change in Qa for a galaxy of some 
arbitrary morphology but averaged over all position angles at this order in Sr/rQ is 

{6qa) ^ hjk J (fr{fo{r))didjdkWa{r). (75) 

This vanishes since {fo{r)) is an even function while didjdkWa{r) is odd, provided we take w{r) to be 
circularly symmetric at least, which we will assume is the case. The effective net fractional error in the 
KSB approximation is therefore on the order of {Sr/ro)'^, or typically on the order of the induced stellar 
ellipticity. 



4.2. General PSF Anisotropy Correction 

In principle, one can develop a higher order correction scheme within the context of this model, but 
this does not seem to be particularly promising; it is not clear that the result will be sufficiently robust and 
accurate for even for ground based observations — in very good seeing conditions the PSF anisotropy from 
aberrations becomes large, and the perturbative approach will break down. Moreover, for observations with 
telescopes in space, the model of the PSF as a convolution of a perfect circular PSF with a kernel, small or 
otherwise, is wholly unfounded. The instantaneous OTF is 

.g(fc) = J (fr A{r)A{r + kD\/2n) exp{i[ip{r) - (p{r + kD\/2n)]) (76) 

where A is the real transmission function of the telescope input pupil, D is the focal length, and (p{r) is the 
phase error due to mirror aberrations and the atmosphere. It is often said that the general OTF factorizes 
into a set of terms describing the atmosphere; the telescope aperture; and aberrations, and this would seem 
to justify the convolution model discussed above. For atmospheric turbulence, and for aberrations arising 
from random small scale mirror roughness this is correct. This is because the average of the complex 
exponential term in (|7^) depends only on Sr = kD\/2'K and is independent of r, and (|7^) then factorizes 
into two independent terms, and the PSF is then the convolution of two completely independent and 
non-negative functions, but this is not the case in general. To be sure, one can write the combined aperture 
and aberration OTF as a product of some 'perfect' OTF with some other function (the true OTF divided 
by the perfect one) but quite unlike the case for random turbulence and mirror roughness, this function is 
neither independent of the shape of the pupil, nor is it positive; if you compute this function for a telescope 
subject to a low order aberration from figure error, for instance, you will find that the this function is 
strongly oscillating and is just as extended as the true PSF. Consider the situation in WFPC2 observations, 
for example, where the PSF anisotropy has important contributions from both from asymmetry of the 
pupil A{r) and from phase errors, though with the former tending to dominate (though not enormously so) 
for long wavelengths and far off axis with the WFPC2. The aperture function for the lower left corner of 



chip 2 computed from the Tiny-Tim model (Krist 1995) is shown in figure The off-axis pupil function is 



approximately that on-axis minus a disk of some radius ri at some distance Ar off-axis. If we let the radius 
of the primary be rp and define a disk function Drg{r) = O(r/ro — 1) then on computing the electric field 
amplitude a{x) and squaring we find that the on-axis PSF is of course given by the Airy disk: g ~ D^^ 
while the off axis PSF g' is given by g'{x) ~ g{x) — Dr^Dro cos{2TTxAr / DX). For ri ^ tq, Dr^ is relatively 
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slowly varying and close to unity, in which case the extra off-axis obscuration introduces a perturbation 
which is proportional to the product of Dq and a planar wave, or essentially an asymmetric modulation of 
the side lobes of the on-axis PSF (see LH panel of figure ||). 

PUPIL AHD PHASE ERROR UaDULAriDM TRANSFER FUNCTIDN 




Fig. 4. — The panel on left shows the WFPC2 pupil function and phase error from Krist's Tiny Tim model 
for an off-axis point on the focal plane. The phase is given in radians for a wavelength of 800 nm. The panel 
on the right shows the modulation transfer function (the real part of the OTF). 

In general then, and especially for diffraction limited observations, the perturbative convolution model 
cannot be trusted. Luckily, at least in the context of the shear estimators discussed here, a simple solution 
clearly presents itself: Since to compute the polarizability requires that one generate a rather detailed model 
of the PSF, and that we probably want to apply some kind of re-convolution fo to obtain a 

well behaved shear response, we might as well use the opportunity to choose g''{r) in order to re-circularize 
the PSF exactly. For example, from the observed PSF, one can compute the OTF g{k) and then form the 
circularly symmetric function gmin(fc) being the greatest real function which lies everywhere below |g(fc)|. 
The function ^^(r) with transform 

3^(fc)=5min(fc)Vff(fc) (77) 

both guarantees a well defined shear polarizability and gives a perfectly circular resulting total PSF. This 
is illustrated in figure |^. 

In some situations the PSF may have near 180 degree symmetry, in which case a good approximation 
is to re-convolve with a PSF which has been rotated through 90 degrees. Since the PSF is real, the OTF 
must satisfy the symmetry g{k) — g*{—k). In the absence of aberrations, a diffraction limited telescope has 
OTF g = A(B A which is real and so has exact symmetry under rotation by tt, a property which is shared 
the PSF, so RT^g — g. If one convolves with g'^ = i?7r/25 then the resulting total PSF g'' ^ g is symmetric 
under 90 degree rotations (proof: i?7r/2fftot = R-,T/2{R-n/2gg) = R-^gRT^/ig = gR-n/2g = 5tot)- For a galaxy of 



-24- 



SOURCE PSF 



RECIRCULARISING RSR 
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1 1 1 1 1 1 1 1 1 

if 
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1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 

|fi 








Fig. 5. — HST WFPC2 PSF from Tiny Tim and re-circularizing filter computed as described in the text. 



^iven form, the average PSF induced polarization is 

(fc) = / d^r {gtot «) f)Wa = [(fftot (8) /) (8) Wa]r=0 



{2nf 



\otWo 



(78) 



where Wa = Maimw{r)rirm- But (/) is circularly symmetric and tlia is anti-symmetric under 90 degree 
rotations, so {qa) = 0. Other situations where this would work would be groimd based observations but 
with large amplitude telescope oscillations, and in fast guiding. In general, however, aberrations from figure 
errors will introduce both real and imaginary contributions to the OTF and the latter will destroy the exact 
180 degree symmetry of the PSF. This is certainly the case for the WFPC2 PSF at optical wavelengths. 

In the linearized approximate schemes described earlier, the PSF anisotropy is entirely characterized 
by the two coefficients ka- In general, these will vary substantially with position on the image, but this 
can be treated by modeling ka as some smoothly varying function of image position R such as a low order 
polynomial. In the scheme proposed here we need to be able to generate not just these two coefficients, 
but the full two dimensional PSF g{r, R) (being the intensity at distance r from the ccntroid of a star 
at position R). A simple practical approach is to solve for a model g{r,R) ~ Yli 9ii''')h{R)^ \vit\i fi{R) 
being some set of polynomial or other basis functions. A least squares fit for the image valued mode 
coefficients is gi{r) = mjjGj{r) with m/j = Estars (^li^) = Estars /j'(^)5obs(?-)- Armed 
with the image valued coefficients one can then, for example, compute the convolution fg ^ g ® fo as 
fsi"!") — Ylii fi{'''){9i ^ fo), that is we convolve the source image with each of components gi{r) and 
then combine with spatially varying coefficients fi{R)- Non-linear functions of the PSF such as the 
re-circularizing kernel g^ are easily generated by making realizations of the PSF models on say a coarse grid 
of points covering the source image, and from each of these computing the desired function, and then fitting 
the results to a low order polynomial just as for g{r, R). 

The approach described here lends itself very nicely to observations with large mosaic CCD cameras, 
in which misalignment of the chip surfaces and the focal plane coupled with telescope aberrations can give 
rise to a PSF which varies smoothly across each chip, but which changes discontinuously across the chip 
boundaries. This results in a very complicated PSF pattern on the final image obtained by averaging over 
many dithered images. It is a good deal simpler to re-circularize each of the contributing images. We still 
need to generate a spatially varying model of the final g{r), g^{r) in order to compute the polarizability, 
but if we fail to model these exactly, it will only introduce a relatively minor error in the shear-polarization 
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calibration. 



4.3. Noise Bias 

In the absence of noise, the recircularization procedure exactly annuls any PSF anisotropy; the signal 
content of the image is exactly statistically isotropic. The noise component in the images (which in the 
original images is incoherent Poisson noise) will however be correlated with an anisotropic two-point 
function - the peaks and troughs of the noise will appear as ellipses with correlated position angles. For a 
given object, the noise is equally likely to produce a positive or negative fluctuation in the polarization g^, 
and the effects cancel out on average. Due to the nature of faint galaxy counts, however, objects of a given 
observed flux are more likely to be intrinsically fainter objects which have been scattered upward in flux 
than intrinsically brighter objects which have been scattered down, and there will therefore be a tendency 
for faint objects to be aligned like the re-circularizing PSF; i.e. oriented opposite to the PSF. We will 
refer to this as 'noise-bias'. It should noted that an analogous effect is present in the old KSB anisotropy 
correction scheme. In that case the noise is isotropic, but objects are detected as peaks of a smoothed 
image and there is then a tendency for the galaxies close to the threshold for detection to be aligned like 



the PSF. In the context of perturbative schemes as described in §4.1 there are other unevaluated errors in 
the correction which are typically of similar order to the noise bias effect, which goes some way to explain 
why it has been ignored in the past. In the improved PSF anisotropy correction method proposed here it is 
the leading source of error and it behooves us to analyze and correct for it. 

Noise will affect the photometric parameters like in two ways; in addition to a straightforward 
additive noise term, noise will also affect the object detection, and will shift the centroid about which we 
measure the centered second moments, for instance. In fact, the de-centering effect is quite weak. To see 
why, consider the following simple (though actually quite practical and realistic) model for the detection 
and measurement process in which we take the re-circularized image fs and define objects as peaks of the 
field F{r) = w (B fs where w is a some smooth, circularly symmetric weight function (we typically use a 
compensated 'Mexican hat' filter), and the object moments qa are the value of the field qAif) — wa ® fs 
evaluated at the peak location. Let us assume that there exists some object which, in the absence of noise, 
would lie at the origin, so for this object we would have 

F ^F{0) = J (frw{r)fs{r) 

0^dm = S<Pr{d,w)fs{r) (79) 
qA = 9^(0) ^ J (Pr w Air) f sir) 
where the condition that the object be a peak of F is expressed as the vanishing of the dipole-like quantity 
diir) = diF. If we now add noise, and let fs^fs+ fm where fnir) is the noise component of the 
re-circularized image, all of the photometric fields Fir), diir), QAir) will change to F' — F + SF ~ F + w® fn 
etc. The dipole di will no longer vanish at the origin, but will have some value diiQ) = J d^r idi'w)fn, 
and the peak of Fir) will have moved to some position rpk ^ 0. In the vicinity of the peak we have 
diir) ~ -|- (r — rpk)j x djdiir) and hence, to first order in the noise amplitude, the peak location is 

-1 



''i.pk = -idjdi) '^Sdj = 



d'^r didjwfs 



d^r djwfn (80) 



The shift of the peak has no effect on F at first order (because F is stationary) but is does affect the central 
moments qa and we have 

SF^Jd^rwir)f^ir) 
6qA = J d'^r [wAir) - VAjdjw] /„ 
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where 



VAi 



d\ id,WA)fs (82) 



Thus, in general, the first order change in the second moments qa depends not just on tlie noise and the 
kernel wa, but on the form of the underlying noise-free image fs through the term involving vaj- There 
is good reason to think that this form dependence is weak for real objects; the function diWA is an odd 
function, so the extra term vanishes if the galaxy is symmetric under 180 degree rotations. Also, it is hard 
to see why the presence of this term would cause a systematic polarization; the change in the polarization 
5qA associated with the centroid shift is proportional to the vector VAj which, being a function of the 
re-circularised image fs is equally likely to be positive or negative. For each galaxy that, for a given 
realisation of noise, suffers a certain Sqa there is a 180 degree rotated clone which has precisely the same 
weighted flux, polarization etc. (these being even functions) but which suffers a Sqa of opposite sign. The 
kind of effect envisioned in the leading paragraph of this sub-section arises if there is a correlation between 
fluctuations in the polarisation and the flux F, coupled with a gradient of the density of objects n{F). 
The first term in Sqa in ( |8l] ) does, as we shall see, indeed correlate with 5F, but, when we average over 
orientation of the underlying noise- free galaxies, the second does not. Unfortunately, we have not been able 
to rigorously demonstrate that the centroid shift term vanishes. Nonetheless, in what follows we will assume 
that this term is negligible. This should certainly be adequate in order to estimate the magnitude of the 
effect, and is probably sufficiently accurate to give a useful correction, though the results should strictly be 
regarded as only approximate. 

Let us now assume one has measured a set of n photometric parameters pi for a galaxy which, like 
F, (Jo, Qa, are linear functions of the brightness /o and of the form Pi = J d^r Ki{r)fo{r). The effect of 
noise in the images will be to introduce a perturbation in these parameters 5pi. These perturbations have 
zero mean, {5pi) — 0, and since there are typically a very large number of photons (noise plus signal) 
contributing to the parameters pi, the central limit theorem dictates that the 5pi will have a multivariate 
Gaussian distribution: 

n{5p)d''6p - ((2^)"|C|)-i/2 exp(-%C^i<5p,/2)rf",5p (83) 
where the covariance matrix is = (SpiSpj) and is given by 



Cij — 



j d^r K,{Kj®ir.) (84) 



where S.ni'i') = {fn{r')fn{'r' + r)) is the 2-point function of the noise /n(r), and which may be computed as 
described in appendix 

For objects which are detected at a fairly high level of significance ^ 1, the noise will cause a small 
modification of the underlying distribution function: n{p) — )■ n'{p), which we can compute in a perturbative 
manner: 

f dn 1 d'^n 

n'{p) = / n{p ~ 5p)a{5p) = n{p) ~ -7^{5p^) + - „ „ {Sp^Spj) + ■■■ (85) 



dpi 2 dpidp 



Since <j{p) is an even function, all the odd terms in this expansion vanish, and we have 
n'{p) = n{j)) -)- ^Cijd^n/dpidpj at leading order. 

Let us specialize now to the n = 4 dimensional case: pi — F,qQ, and compute the mean polarization 
for galaxies of given F, qo with some weight function W{q — \qa\)- 



, . ^ f(Pq W{q)qan'{F, go, qg) J d^q W{q)qaSn{F,qo,qc,) 

[q.)F,,o jcPqW{q)n'{F,qo,qg) ^ J d^ q W [q^F, qo , q^) ^ ' 
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where 6n{F,qo,qa) — ^Cijd^n/dpidpj. Since qa is an odd function and W{q), n{F,qo,qa) are even, the 
only terms which contribute to the integral in the numerator are those involving a single derivative with 
respect to polarization: d^n/dFdqp and d^n/dqodqp and, on integrating by parts, we have 



-Ud'q nSFSqp)§^ + {SqoSqf,)§^ 



(87) 



! <PqW{q)n{F,q^,q^) 

Using d{W (ci)qa) I dqp — 5ai3W{q) + qaqfidW / dq, integrating over angle and converting from n{F,qo,qa) to 
n'{F, qo,q) = 2Trqn{F, qa,qa) and dropping the prime we have 

-i/dq ({5F5q^)^ + {5q^5q^)§£\ {W {q) + \qW' {q)) 

^'-^ - ^ SdqnWiq) 

and letting W{q) become a delta-function to isolate a single value of q we have 

n{qc.) = -\ [{5F5q^) ^ + (SqoSq^)^^ (89) 

Thus, as anticipated, there is a net induced polarization if there is a non-zero correlation between the 
polarization fluctuation dqa and F and/or qo, and also significant gradients of the distribution function 
with respect to F ot qo. 

The effect on the re-scaled polarization of the first term in (g9|) )/F = 

— {^{6FSqa)/F'^)d\nn/d/lnF and is second order in the inverse significance: (q^) cx where 
^2 _ p"^ I l^^p"^"^ and rapidly becomes small for well observed objects. Thus it should be possible to set 
a sensible limit on the significance of object detection, and, if necessary, use ( |89| ) together with (jsj) for 
(SFSqa) etc. to correct for this. The error in this linearized approximation is of fourth order in the inverse 
significance. As a simple illustrative example consider a Gaussian galaxy of scale re, a Gaussian window 
function w with scale r^,, and a Gaussian ellipsoid PSF g = exp(— (x^/r^ -t- ?/^/r^)/2) with ra = rg{l + e/2), 
fb — fgi^" c/2). A suitable recircularising kernel is then — exp(— (a;^/r^ + ?/^/r^)/2) and the normalised 
two point function of the noise is then 

TTran 

The expectation averages are, to first order in e, (SF^) — 2(r^ -I- r^); (SqiSF) — r^{r1^ ~ ^6)/('"^ + """I) = 
2er^rg/(r^ -|- r^), and so taking only the first term in ( p9| ) for simphcity 

using 91nn/91nF ~ 1 as observed. A shear of strength 7 applied to the Gaussian galaxy produces a 
polarisation 

{K)^^irlrl/(rl + rl+rl) (92) 
therefore a PSF anisotropy of strength e is equivalent to an effective shear 



' 0..2 „2 /„2 , 2\2 V^"^/ 

or, for = tg = V^fg say, 7 = 25e/144z^^, which for 1/ = 6 and PSF asymmetry e = 0.3, which is a 
reasonable value for off-axis points on the CFHT in good seeing, gives a shear of around 1.4%. This is a 
sizable effect, and should therefore be corrected for. 
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5. Discussion 

We have considered the problem of how to estimate weak gravitational shear from observations which 
have been degraded by atmospheric and/or instrumental effects. Previous analyses of this problem have 
made simplifying assumptions which render the results inaccurate. A major result of the paper is the 
finite resolution shear operator ( ^ ) which gives the response of an observed image to a gravitational shear 
applied before smearing with the PSF. This result can be used to properly calibrate the effect of any shear 
estimator, and is valid for arbitrary PSF, be it turbulence or diffraction limited. We then focused on the 
application to weighted moment shear estimators. We have computed the response of individual objects 
to a shear in j|3.l|, and the response of the population of background galaxies with given photometric 



properties in § |3.2| , and from this we have devised an optimal weighting scheme §3.3. In the last section we 
have considered the correction for PSF anisotropy. While there are still some approximations in the present 
analysis, we feel that they place the techniques of shear measurement on a much firmer footing than before. 
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A. Properties of Optical Point Spread Functions 

Here we shall briefly review and derive some properties of telescope point spread functions which 
are used above. For more detailed background the reader should consult (Roddier 1981; Beckers 1993| ) 



and references therein. We will highlight the wavelength or color dependence of the various sources of 
PSF anisotropy, which may be crucially important for weak lensing searches for large-scale structure and 
galaxy-galaxy lensing. 



According to elementary diffraction theory ( Born fc Wolf 1964 ) the complex electromagnetic field 



amplitude a{x) due to a distant source at position Xphys on the focal plane (we will suppress polarization 
subscripts for clarity) is given as an integral over the input pupil 

a(a:phys) = / d'r A{r)C{r)e'^'^--y''-/''^ (Al) 



where A{r) is the 'pupil function' describing the aperture transmission, C(r) is the complex electric field 
amplitude of the incoming wave, A is the wavelength of the radiation and L is the focal length. The field 
amplitude C(r) will incorporate any random amplitude and phase variations of the incoming wavefronts 
due to atmospheric turbulence, whereas constant wavefront distortions due to aberrations in the optical 
elements of the telescope are incorporated as a complex factor in A. Thus a{x) is the Fourier transform 
of AC, evaluated at wave-number k = iirix/DX. The PSF g is the square of the field amplitude and in 
rescaled coordinates x = 2nXpi^ys/ LX, is 

gix) = \a{x)\' = J ff(z) (A2) 

where the OTF is 

g{z)^ J d^rC{r)C*{r + z)A{r)A*{r + z). (A3) 

For very short ground-based observations the atmospheric rippling is frozen and the PSF consists of 
speckles. For long exposures we are taking the time average of the OTF and can replace C{r)C*{r + z) 
by its time average {C{r)C*{r + z)) = Cc(-z)- This is independent of r, so the OTF factorizes into two 
independent functions 

g{z)^^c{z) I d^rA{r)A*{r + z) (A4) 



and the same is true for random small scale amplitude or phase fluctuations introduced by e.g. random fine 
scale mirror roughness. 



A.l. Atmospheric Turbulence 



Ground based observations on large telescopes are usually limited by atmospheric seeing arising from 
inhomogeneous random turbulence, and it is a good approximation to ignore the finite size of the entrance 



aperture and set the factor involving A in (A4) to unity. In the 'near field' immediately behind the turbulent 
layer, the effect on the incoming wave is a pure phase shift C(r) = e"^^''^ where (p = 27rid{r)/\ and d{r) is 
the vertical displacement of the wavefront due to turbulence. The displacement d is nearly independent of 
wavelength, so ip{r) oc 1/A. At greater depths this phase shift evolves into a combination of amplitude and 
phase variations, but the 2-point function {C(r)C*{r + z)) remains invariant (Fried 1966; Roddier 1981)and 
the 'natural seeing' OTF is 



J{ip(r)-ip(r+z))\ 



(A5) 
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For steady turbulence and long integrations the central limit theorem guarantees that the phase error 
■0 = ip(^r) — ip{r + z) will have a Gaussian probability distribution p{ijj) = {271 {tp"^))^^^^ exp{—tp'^/2{tp'^)) and 
so the time average of the complex exponential is 



(e'^) ^Jd^ p{^)e'^ = exp(-(V;2)/2) = exp(-5^(r)/2) (A6) 
where the 'phase structure function' is 

5^(Ar)^((^i~^2n. (A7) 

There are strong theoretical jTatarski 196l|) and empirical reasons to believe that on scales much less than 
some 'outer scale' the turbulence will have the Kolmogorov n = —11/3 spectrum, for which Sdir) oc r^/'^. 
The structure function for the phase is conventionally written as S^{r) = 6.88{r/ro)^^^ where ro is the 
'Fried length' ( [Fried 196^ ) being on the order of tens of cm for typical observing conditions (an ro of 20cm 



gives a FWHM — 0".5 at A = 550 nm). The rms phase difference rises with separation as r^/^ in the 
'inertial range' delimited at the upper end by the outer scale, set by the width of the mixing layer, which 



recent estimates ( Avila et al. 1997 ; Martin et al. 1998 ) find to be around 10 — 20m and much larger than 



rg. The on-axis OTF computed from stars in deep CFHT imaging agrees quite well with the theoretical 
expectation. The inertial range is limited at the low end by diffusion, but at scales much smaller than rg, 
so little error is incurred in ignoring this; in real telescopes mirror roughness and other effects modify the 
OTF at small scales. These effects dominate the PSF at very large radii, but are unimportant for weak 
lensing observations. 

The optical transfer function is g{k) = cxp{—S^{kDX/2TT) /2) and is real and positive. The PSF is the 
transform of exp(— 6.88(2:/ro)^ '^) with width which scales as i?FWHM oc A~^/^ which again is found to apply 
quite well in practice. This very weak dependence on wavelength is a blessing in weak lensing since one 
uses stars to measure the PSF for the galaxies, yet the stars and the galaxies may have different colors. 
The atmospheric PSF is expected to be isotropic. At large angles the PSF has profile g cx x"^^^"^ so the 
unweighted second moment of the PSF is not well defined. For Kolmogorov turbulence the log of the OTF 
is just proportional to k^^^ and is well defined for all k. 



A. 2. Fast Guiding 



According to the Kolmogorov law, the rms wave-front tilt, averaged over scale r varies as r which 
suggests that even for telescopes with reasonably large D/ro there may be useful gain in image quality from 



fast guiding, and experience with HRCAM on CFHT (McClure et al. 1989) would seem to support this, 
though part of the dramatic improvement is likely due to inadequacy of the existing slow guiding system. 
A technological advance which may have implications for weak lensing is the advent of on-chip fast guiding 



( Tonry, Burke, fc Schcchtcr 1997 ) with OTCCD chips. With a mosaic camera composed of such devices it 
should be possible to obtain partial image compensation over a large angular scale, with one or more guide 
stars for each 'isokinetic' patch. 



The theoretical fast guiding PSF was first explored by Fried (1966 ) who argued that the OTF should 
take the form of the natural seeing or uncompensated OTF times an 'inverse Gaussian' exp(-|-Q!fc^), with 
scale factor a given in terms of D and tq. This is a physically reasonable picture, since it implies that the 
natural PSF is the convolution of the corrected PSF with a Gaussian to describe the distribution of tilt, but 



is only an approximate result. Modified forms of the 'Fried approximation' have been explored by Young 
(1974) and Jenkins (1998| ), and the fast guiding OTF has been simulated by Christou (1991). 
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It can be shown that in the near-field hmit the exact fast-guiding OTF is given by 

g(z) = / (frA{r)A{r + z)exp{^{tl;{r,zf)/2) 



where 



{-^{r, zf) = S{z) + z,[{W, ® S)r - {W, <E) S), 



1 



(fr' Wi{r'){Wj S)r 



(A8) 



(A9) 



and where Wi = di A^. Examples are shown in figure These plots show that the impact of fast guiding 
on the atmospheric PSF for large telescopes will be rather modest, at least if current, rather low, estimates 
of the outer scale are correct. Fast guiding may, however, yield dramatic improvements for small Im 
diameter) telescopes. How well this would work depends largely on the altitude of the turbulent layer. The 
isokinetic patch size (over which stars move coherently) is ~ Z?//i so for h ~ 10km and D = 1, say, this is on 
the order of 20", the motion needs to be sampled at a rate ^ v/D reflecting the relatively high wind speed 
at high altitude, and it may then be difficult to find bright enough guide stars. There are strong indications 



( Christian fc Racine 1985 ; Tonry, Burke, fc Schechter 1997 ; McClure et al. 1991 ) that centroid motions 
are coherent over much larger angular scales than this, indicating that much of the image degradation 
arises from low-altitude turbulence, and this greatly improves the outlook as one can determine the local 
motion by averaging a number of stars, and one can afford to sample at a lower rate. A collection of small 
telescopes equipped with wide angle OTCCD cameras could be a formidable instrument for weak lensing or 
other projects requiring high resolution imaging over wide fields. Fast guiding, while offerering important 
resolution gains, will also present its own challenges since one expects the PSF to become sytematically 



anisotropic depending on location with respect to the guide stars for the reasons described by McClure 



et al. (1991). Also, fast guiding does not cure PSF anisotropics from telescope aberrations. 
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Fig. 6. — Optical transfer functions and corresponding PSFs given by equations (A8, A9) for a range of 
telescope aperture diameters: 1.0m solid; 1.5m dash; 2.2m dot-dash; 3.6m dotted; oo dot-dot-dot-dash. A 
von Karman turbulence spectrum with outer scale of 20m and Fried length tq = 0.2m were assumed, and 
the telescope was assumed to be operating at a wavelength A — 550nm. 



A. 3. Atmospheric Dispersion 



Variation of the refractive index of the atmosphere with wavelength will cause images to be dispersed 
into spectra, and the consequent elongation of images was estimated by KSB, who pointed out that this 
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could potentially cause problems when using PSF's of stars to correct the shapes of galaxies, since the 
elongation depends on the spectrum of the object; an object with a line-like spectrum, or one with a 
sharp-edged absorption band which falls within the bandpass, will of course be less dispersed than a 
continuum object. Here we will make this more quantitative. 

Weak-shear measurements are usually made on images composed from multiple exposures taken at 
a range of air masses and which are co-registered using centroids of foreground stars. The astrometric 
solution obtained using stars of a variety of spectra will be a compromise. Assuming, as is usually the case, 
that one has many more stars at ones disposal than the number of coefficients of the transformation one 
is trying to solve for (typically a low order polynomial), then the solution obtained by minimising least 
squared residuals will be one which correctly maps surface brightness at some wavelength Aq determined 
by the average properties of the foreground stars, but with photons at other wavelengths displaced by an 
amount proportional to A — Aq . The sum of a set of N such images is 

N 

/tot(r,A)=^/(r + a,(A-Ao),A) (AlO) 
1=1 

where f{r,X) is the image seen at zenith. The 2-vector an has \ai\ = (0o/Ao)(91n8/91nA)tanz7, with zi 
the zenith distance, and ai is directed towards the horizon. Tables of the refraction angle 6(A) are given 



by Allen (1973 ). A star with SED f\, and therefore with photon distribution function S{X)dX oc Xf\dX, 
gives, in a photon counting system, a response at zenith /(r. A) = S{r)S{X), and the summed image of such 
a star is then 



N 



/tot(r) = f dA/tot(r, A) = ^ / dXS{r - aj{X - Ao))5(A) (All) 

1=1 

The centroid of a star in the image /tot (t) is an average of the centroids that would be obtained from each 
of the contributing images. Consider the /'th image, and rotate the coordinate system so that a/ lies along 
the X-axis. The centroid for this component is r/ = (x/, 0) with 

_ ^ JdxJdyJdX xSjx - aj(A - Xo))S{y)S{X) ^ aj J dX {X- Xo)S{X) ^ (T _ \ \ 

JdxJdyJdX6ix-aiiX-Xo))6iy)SiX) J dX S{X) "'^ °' ^ ' 

from which it follows that the Aq which minimises (x/), being the average over the stars used for registration 
of the squared displacement, is Aq = (A) = — ^ A, i.e. simply the average over the registration stars of 
their mean wavelength. The centroid for the summed image is 

n = (A- Ao)(a.)/ (A13) 

where (a;)/ = j^'^ aii- Similarly, the central second moment is 

/dV {r,~ri){rj~rj)ftot{r) 



= \,2 V , ^ = - Ao)^(«.a,), - (A - X,Y{a.)i{a,)i (A14) 

Note that since, for any sensible zenith angle, the width of the spectrum is tiny compared to that of the 
PSF, this unweighted second moment fully characterises the effect of dispersion. 

Consider, for illustration, the case of an equatorial field observed with a telescope near the equator. 



and for a range of zenith angle \z\ < Zmax- In this case ^ a/ = 0, so the second term in (A14) vanishes 
(this is not a particularly unusual special case; for a number of observations spread over a range of zenith 
angles, one would generally expect the second term here to become relatively small). When we solve for the 
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PSF from the shapes of foreground stars of various types we are effectively averaging over a mix of SED's 
appropriate for a low redshift spiral galaxy. The average PSF moment obtained from the stars can, for the 
equatorial case, then be written as 

(p..) = ((A^>-(A)2)(a^)7 (A15) 

where by (a^)/ = (9o/Ao)^(<91n6/91n A)^;^ tan^ 2/, and the SED dependent coefficient of {0^)1 is the 
mean of the dispersions of the stars plus the dispersion of the means of the stars. For the equatorial case 
this is the same as the second moment for a single point like object with an SED like the average SED 
of the registration stars. This is not strictly true in general, but the difference is probably a minor one, 
and it then follows that if the faint galaxies have SED's such that their (A — Aq)^ differs systematically 
from that for a low redshift spiral then the PSF correction will be systematically in error. Fi gure |^ 



quantifies the importance of this effect, by using redshifted galaxy SED's of various types from Coleman 



Wu, fc Wecdman (1980| ). The panel on the left shows the displacement of the centroid of galaxies of 



the various indicated types as their spectra are red-shifted, and the right hand panel shows the variation 
of the second moment of the PSF. We find that the Tband second moment is very stable, with peak 
fluctuations 5pim ~ SOOmas^ for these parameters, and with little systematic difference from a low spiral 
redshift galaxy if we integrate over a range of redshifts. For illustration, if we take the systematic change 
in polarisation to be say Spim ^ lOOmas^ and a Gaussian profile galaxy, this corresponds to a shear of 
7 ~ 1.38Jpi™/FWHM^ ~ 6 X 10-''(0".5/FWHM)2. In the V-band the polarisation fluctuations are larger 
by a factor 2-3, but even then, and even for marginally resolved objects in excellent seeing the effect is 
at or below the sensitivity for current and near future surveys. The effect scales as (tan^ z) however, so 
observations at z ^ 1 should be avoided. This would also become more of an issue with fast guiding, where 
we can hope to resolve much smaller objects. The effect is also dependent on the details of the system 
response function; the CFH12K system having a particularly broad I-band response which exacerbates the 
effect. 

In the foregoing we assumed that the foreground stars used for registration were numerous and well 
mixed. With a finite number of stars there will be some additional color dependent error in the registration 
which, if uncorrected, would give rise to artificial field distortion. However, we have found from simulations 
using the USNOA catalog as the astrometric reference system that this introduces artificial shear at around 
the 10^* level, which is negligibly small for present and near future surveys. 



A. 4. Aberrations 

Aberrations of the optical elements of the telescope can be a significant contribution to the anisotropy 
of the PSF. These can be analyzed in the same manner as the wavefront deformation due to the atmosphere, 
but with a couple of distinctive features: First, for low-order 'classical' aberrations where the phase error 
varies smoothly, and for ground based observing conditions, if an aberration is an important factor then 
it's contribution to the OTF will be nearly achromatic. This is because if there is a smooth variation of the 
wavefront error (rather than a Gaussian random field with power at all scales as in atmospheric turbulence) 
amounting to iV 3> 1 wavelengths, then the PSF will be very well approximated by its geometric optics 
limit with shape defined by the pattern of caustics (though for a narrow band filter, the PSF would actually 
be found on close examination to be composed of a set of speckle sized patches concentrated along the lines 



where the classical caustics form (Berry & Upstill 1980)). The wavefront deformation can be measured 



directly from out of focus images ( Roddier fc Roddier 1993| ) so this contribution to the PSF can be directly 
predicted. 
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Fig. 7. — Atmospheric dispersion. Left hand panels show the angular displacement due to changes in the 
SED for galaxies at a range of redshifts. A zenith angle of 45 degrees was assumed, and upper and lower 
panels show the deflection for the I-band and V-band system response functions for the CFH12K camera. 
More concretely, the quantity plotted is the mean wavelength A — Aq = J dX {X — Xo)XfxR{X)/ J dX Xf\R{X) 
where fx is the SED, R{X) the system response, Aq was taken to be 820nm, and A — Aq was converted 
to angle as described in the text, and assuming a 4000m observatory altitude. The V-band plot has been 
limited to z < 1.5 since the CWW SED's are limited to proper wavelength A > 200nm and redshift past the 
y-band at around z = 1.5. The plots shows that the shifts due to SED variation are typically on the order of 
20mas in the I-band, but somewhat larger in V. The elliptical/SO SED shows a somewhat greater excursion 
at high redshift, but this it caused by the SED red-shifting right out of the filter band, so one would not 
expect to find many such objects in flux limited samples. The panels on the right show how the width of the 
spectrum, and hence the PSF polarization, varies with redshift and galaxy type. Here the quantity plotted 
is (A — Ao)^ = J dX {X — Xo)^XfxR{X)/ J dX XfxR{X), again converted to angle, here assuming (tan^ z) = 1. 
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A. 5. Diffraction Limited Seeing 

Truly diffraction limited seeing arises when the rms phase error across the aperture (due to the 
atmosphere and/or aberration of the optical elements of the telescope) is much less than unity. In this case 
the optical transfer function g{k) is, to a good approximation, just the auto-correlation of the aperture 
A(B A a.t lag Ar = kDX/2Tr, and must therefore vanish for spatial frequencies k > 27r/(/A), where / is 
the ratio of the aperture diameter to the focal length. The log of the OTF becomes ill defined as one 



approaches the diffraction limit. From ( A4) we see that this cut-off is also present in the case of turbulence 
dominated seeing, but it occurs at a high frequency where the the optical transfer function has already 
become exponentially small due to atmospheric effects, and has little impact. 

In the absence of aberrations, the OTF for diffraction limited seeing is real and non-negative. The 
OTF is symmetric under rotations of 180 degrees, and so any quadrupole anisotropy of the PSF can be 
anulled simply by re-convolving one's image with a 90 degree rotated PSF. For diffraction limited seeing 
the size of the PSF scales as the inverse of the wavelength, a fact which can be incorporated in empirical or 



theoretical ( Krist 1995 ) modeling of the PSF 



In the HST WFPC2, figure errors are not negligible. The imaginary part of the OTF is excited to 
the degree that re-convolution with a 90-degree PSF still leaves non-negligible PSF anisotropy. The phase 
errors are not large — the telescope is nearly diffraction limited — so this means that the wavelength 
dependence could be quite complicated. The systematic error arising from differences between faint galaxy 
and foreground star SED's can be estimated much as we did for atmospheric dispersion. 



A. 6. Guiding Errors, Pixellisation, and Detector Effects 

So far we have considered the continuous distribution of intensity on the focal plane /ol?")- In real 
detectors we sample the image with a grid of pixels. The response of a pixel in not uniform and has 
been directly measured for front-illuminated EEV devices by scanning a small spot of light across the 



CCD (lorden, Deltorn, & Oatcs 1993, Jordcn, Deltorn, & Oates 1994). They found very little 'leakage' 



of electrons across pixel boundaries, but according to the Krist (1995) this is a substantial effect for the 
WFPC2 instrument on HST; a photon landing in one pixel has a non-negligible probability of being detected 
by a neighboring pixel instead. Such effects are expected, and found, to depend on wavelength. The value 
of a pixel is a sample of the convolution of the sky surface brightness with the 'pixel function' p(r). 

In many real systems the pixel spacing d is not much less than the instrumental resolution and images 
from single exposures are quite badly under-sampled. For this and other reasons, images are typically 
constructed from a series of exposures with either systematically (in the case of HST) or randomly (for 
terrestrial observations) staggered positions on the sky. Each image gives a 2-dimensional grid of samples of 
p® foi and a piecewise continuous function // can be constructed by shifting the grid of delta- functions into 
an absolute astrometric coordinate system and convolving with some interpolation function pinterp- One 
can incorporate the effect of guiding errors on a single exposure as a convolution with the pixel function. If 
we average a set of N such images the result is 

F{r) = ^ ^[(/o ® p).CaA ® Pintcrp (A16) 
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where ca(t") is a 2-dimensional comb function with spatial offset (in units of the pixel spacing) A: 

oo 

CA(r)= Yl Hr-{i + A)) (A17) 

ix,iy=—oo 

The form of pinterp depends on the type of interpolation used. For 'nearest pixel' interpolation pinterp is just 
a uniform box of side d, but if one linearly interpolates between the pixel samples, for example, then Pinterp 
will be a more extended, but again readily computable, function. 

The transform of F{r) is 

oo ^ 

-F(fc) = Pinterp ^ (/oP)fc-27rm/d^ exp(27rim • A/) (A18) 

■mx,Tny = — oo I 

The rrix = niy = term in the first sum here is just the ideal image /o convolved with p and with pinterp, 
while the higher order terms represent aliasing. The transform F{k) being Pinterp times the superposition 
of a grid of images of p fo with spacing 27r/(i. Since p is a fairly compact function the dominant aliasing 
comes from the low-order images m = ±1. Aliasing is most severe for a single exposure since the low-order 
aliased images contribute with unit weight. If we average N randomly shifted images then the strength of 
the m ^ terms is reduced by a factor ^ 1/y/N, aliasing is greatly reduced, and to a good approximation 
the field F{r) is simply the convolution of /□ with pCSpintorp- With systematically staggered images, as is 
possible with HST and potentially with fast on-chip guiding, one can do even better; with a uniform M x M 
grid of offsets covering the unit pixel, the nearest, and therefore most problematic, modes rn^, my = ±1 are 
then zero and the modes remain small until we get to a multiple of 2M times the Nyquist frequency. This 
assumes that the transformation from detector to sky coordinates is determined and applied accurately. If 
we make finite errors 6A in registration, the resulting image will be the convolution of the ideal PSF with 
a highly compact cluster of delta- functions, and the optical transfer function g will be the product of the 
atmospheric, telescope transfer functions with the Fourier transform of this pattern. In practice, one can 
typically register images to a small fraction of a pixel (say ^ 0.05 pixels), and the effect of inaccuracy at 
this level will have negligible effect on the final PSF. 

Noise in the images, assumed to be incoherent Poisson noise in the source images, can be analyzed in a 
similar manner and we find that the two-point function of the noise is just the convolution of Pinterp with 
itself, and the two-point function of the noise in re-circularized images can be obtained by convolving the 
raw noise ACF with twice. 

B. Simulated Data 

To test the procedures described here we have generated simulated mock data and then analysed these. 
The simulations were made to match as closely as possible observations of ~ 3hr integration on the CFHT 
with 0".6 seeing. 

We first generated a set of 200 mock catalogue of galaxies each corresponding to a patch of sky of size 
2'. 56 on a side. Galaxies were drawn from a Schechter style luminosity function laid down in a Poissonian 
manner in an Einstein de Sitter cosmology. Images with pixel scale 0".075 were then generated by realising 
the galaxies as exponential disks with random orientations and scale lengths corresponding to fixed rest 
frame central surface brightness. The galaxies were modelled as optically thick, since the optically thin 



- 38 - 



model looks unrealistic as it has too many very bright edge-on systems as compared to real images. A 
number of point-like stars were added to the images, which were then sheared with 7 = 0.1, convolved 
with a Kolmogorov turbulence PSF with 0".6 FWHM, and then rebinned to 0".15 pixel scale. When the 
real data are analysed they are interpolated from the original 0".2 pixel scale to the final 0".15 scale with 
bi-linear interpolation. This results in a further convolution of the signal and the noise in the real images, 
but with slightly different smoothing kernels. These kernels were computed by modelling the image shifts 
as a uniform distribution within the final pixel size; the noise-free mock images were convolved with the 
appropriate kernel and then Gaussian white-noise images were generated to model the sky noise, and were 
convolved with the appropriate kernel and then added to the images. A sample image is shown in figure 
^, and the corresponding size-magnitude diagram is shown in figure |9[ from which it is apparent that the 
simulated objects have properties very similar to those detected in the real data. 

SIMULATED IMAGE REAL CFHT IMAGE 




Fig. 8. — Simulated and real image sections 2'. 56 on a side. 



These data were analysed exactly like the real data. That is, the objects were detected as peaks of a 
smoothed image. The stars were extracted and their shapes fit in the manner described to obtain the PSF. 
A smoothed image fs=9®fo was generated and from this the polarisation was computed using (45), 
and the polarisability for each object was computed using (48). 
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Fig. 9. — Catalogs from image sections shown in figure 



